Bayesian model comparison

1 Setup

Code

rm(list=ls())
libs<-c("rstan", "gdata", "bayesplot", "stringr", "dplyr", "ggplot2", "loo", "hBayesDM", "tidyr")
sapply(libs, require, character.only=TRUE)

Loading required package: rstan

Loading required package: StanHeaders

rstan version 2.26.22 (Stan version 2.26.1)

For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
change `threads_per_chain` option:
rstan_options(threads_per_chain = 1)

Loading required package: gdata

Attaching package: 'gdata'

The following object is masked from 'package:stats':

    nobs

The following object is masked from 'package:utils':

    object.size

The following object is masked from 'package:base':

    startsWith

Loading required package: bayesplot

This is bayesplot version 1.11.1

- Online documentation and vignettes at mc-stan.org/bayesplot

- bayesplot theme set to bayesplot::theme_default()

   * Does _not_ affect other ggplot2 plots

   * See ?bayesplot_theme_set for details on theme setting

Loading required package: stringr

Loading required package: dplyr

Attaching package: 'dplyr'

The following objects are masked from 'package:gdata':

    combine, first, last, starts_with

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

Loading required package: ggplot2

Loading required package: loo

This is loo version 2.7.0

- Online documentation and vignettes at mc-stan.org/loo

- As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session. 

Attaching package: 'loo'

The following object is masked from 'package:rstan':

    loo

Loading required package: hBayesDM

Loading required package: Rcpp

This is hBayesDM version 1.2.1

Attaching package: 'hBayesDM'

The following object is masked from 'package:bayesplot':

    rhat

Loading required package: tidyr

Attaching package: 'tidyr'

The following object is masked from 'package:gdata':

    starts_with

The following object is masked from 'package:rstan':

    extract

    rstan     gdata bayesplot   stringr     dplyr   ggplot2       loo  hBayesDM 
     TRUE      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE 
    tidyr 
     TRUE 

Code

datapath <- '/fast/work/groups/ag_schlagenhauf/B01_FP1_WP2/WP2_ILT_CODE/02_Behav_and_Comp_Modeling'
out_path <- '/fast/work/groups/ag_schlagenhauf/B01_FP1_WP2/WP2_ILT_CODE/02_Behav_and_Comp_Modeling/Output'
behavpath <- '/fast/work/groups/ag_schlagenhauf/B01_FP1_WP2/ILT_DATA'

cbPalette <- c( "#0072B2", "#D55E00", "#009E73", "#CC79A7",
                "#F0E442", "#56B4E9", "#999999", "#E69F00")

loadRData <- function(fileName){
#loads an RData file, and returns it
    load(fileName)
    get(ls()[ls() != "fileName"])
}

2 Step 1: Comparison of random-intercept models

2.1 Load Stan fit files

Code

# stan input non-hierarchical
stan_data <- loadRData(file.path(behavpath,"Input/stan_data_bandit2arm_delta_PH_withC_n58.RData"))

# stan input hierarchical
stan_data_hierarchical <- loadRData(file.path(behavpath,"Input/stan_data_bandit2arm_delta_PH_withC_hierarchical_group_n58.RData"))

fit_file_main <- 'fit_n58_2024-04-12_bandit2arm_delta_main_estimation1_delta0.9_stepsize0.5.rds'
fit_file_main_DU <- 'fit_n58_2024-04-12_bandit2arm_delta_main_DU_estimation1_delta0.99_stepsize0.1.rds'
fit_file_PH_withC <- 'fit_n58_2024-05-07_bandit2arm_delta_PH_withC_DU_estimation1_delta0.999_stepsize0.1_treedepth12.rds'
#fit_file_PH_withC_DU <- 'fit_n58_2024-04-12_bandit2arm_delta_PH_withC_DU_estimation1_delta0.99_stepsize0.1.rds'

fit_main <- readRDS(file.path(out_path, fit_file_main))
fit_main_DU <- readRDS(file.path(out_path, fit_file_main_DU))
fit_PH_withC <- readRDS(file.path(out_path, fit_file_PH_withC))
#fit_PH_withC_DU <- readRDS(file.path(out_path, fit_file_PH_withC_DU))

2.2 Extract log likelihood

Code

# extract log likelihood for each choice
log_likelihood_main <- extract_log_lik(fit_main, parameter_name = "log_lik", merge_chains = T)
log_likelihood_main_DU <- extract_log_lik(fit_main_DU, parameter_name = "log_lik", merge_chains = T)
log_likelihood_PH_withC <- extract_log_lik(fit_PH_withC, parameter_name = "log_lik", merge_chains = T)
#log_likelihood_PH_withC_DU <- extract_log_lik(fit_PH_withC_DU, parameter_name = "log_lik", merge_chains = T)

# create logical vector coding if a column should be included in log_likelihood
x <- logical(stan_data$T) # create vector containing total number of trials * FALSE
a <- c(stan_data$N+1:(stan_data$N*49)) 
x[a] <- TRUE # set logical value to TRUE if index in a

# exclude 1st trial per subject as likelihood is identical and hinders loo estimation (causes Inf pareto k values)
log_likelihood_main <- log_likelihood_main[,x] 
log_likelihood_main_DU <- log_likelihood_main_DU[,x]
log_likelihood_PH_withC <- log_likelihood_PH_withC[,x] 
  
# exclude missing trials
log_likelihood_main <- log_likelihood_main[,log_likelihood_main[1,]!=-999]
log_likelihood_main_DU <- log_likelihood_main_DU[,log_likelihood_main_DU[1,]!=-999]
log_likelihood_PH_withC <- log_likelihood_PH_withC[,log_likelihood_PH_withC[1,]!=-999]
#log_likelihood_PH_withC_DU <- log_likelihood_PH_withC_DU[,log_likelihood_PH_withC_DU[1,]!=-999]

2.3 Fit per model

Code

loo_main <- loo(log_likelihood_main)

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

Code

print(loo_main)

Computed from 36000 by 5583 log-likelihood matrix.

         Estimate   SE
elpd_loo  -1497.2 40.6
p_loo       134.2 12.7
looic      2994.5 81.3
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume independent draws (r_eff=1).

Pareto k diagnostic values:
                         Count Pct.    Min. ESS
(-Inf, 0.7]   (good)     4935  88.4%   228     
   (0.7, 1]   (bad)       425   7.6%   <NA>    
   (1, Inf)   (very bad)  223   4.0%   <NA>    
See help('pareto-k-diagnostic') for details.

Code

plot(loo_main)

Warning in plot.psis_loo(loo_main): 0.29% of Pareto k estimates are Inf/NA/NaN
and not plotted.

Code

if (any(pareto_k_values(loo_main) > 0.7)) {
  loo_main_clean <- loo(log_likelihood_main, k_threshold = 0.7)
  print(loo_main_clean)
  plot(loo_main_clean)
}

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Can't fit generalized Pareto distribution because all tail values are
the same.

Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

Computed from 36000 by 5583 log-likelihood matrix.

         Estimate   SE
elpd_loo  -1497.2 40.6
p_loo       134.2 12.7
looic      2994.5 81.3
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume independent draws (r_eff=1).

Pareto k diagnostic values:
                         Count Pct.    Min. ESS
(-Inf, 0.7]   (good)     4935  88.4%   228     
   (0.7, 1]   (bad)       425   7.6%   <NA>    
   (1, Inf)   (very bad)  223   4.0%   <NA>    
See help('pareto-k-diagnostic') for details.

Warning in plot.psis_loo(loo_main_clean): 0.29% of Pareto k estimates are
Inf/NA/NaN and not plotted.

Code

loo_main_DU <- loo(log_likelihood_main_DU)

Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

Code

print(loo_main_DU)

Computed from 36000 by 5583 log-likelihood matrix.

         Estimate   SE
elpd_loo  -1515.1 42.0
p_loo       122.9 13.2
looic      3030.2 83.9
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume independent draws (r_eff=1).

Pareto k diagnostic values:
                         Count Pct.    Min. ESS
(-Inf, 0.7]   (good)     4847  86.8%   899     
   (0.7, 1]   (bad)       340   6.1%   <NA>    
   (1, Inf)   (very bad)  396   7.1%   <NA>    
See help('pareto-k-diagnostic') for details.

Code

plot(loo_main_DU)

if (any(pareto_k_values(loo_main_DU) > 0.7)) {
  loo_main_DU_clean <- loo(log_likelihood_main_DU, k_threshold = 0.7)
  print(loo_main_DU_clean)
  plot(loo_main_DU_clean)
}

Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

Computed from 36000 by 5583 log-likelihood matrix.

         Estimate   SE
elpd_loo  -1515.1 42.0
p_loo       122.9 13.2
looic      3030.2 83.9
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume independent draws (r_eff=1).

Pareto k diagnostic values:
                         Count Pct.    Min. ESS
(-Inf, 0.7]   (good)     4847  86.8%   899     
   (0.7, 1]   (bad)       340   6.1%   <NA>    
   (1, Inf)   (very bad)  396   7.1%   <NA>    
See help('pareto-k-diagnostic') for details.

Code

loo_PH_withC <- loo(log_likelihood_PH_withC)

Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

Code

print(loo_PH_withC)

Computed from 36000 by 5583 log-likelihood matrix.

         Estimate   SE
elpd_loo  -1377.0 38.5
p_loo       152.2 10.9
looic      2754.0 77.0
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume independent draws (r_eff=1).

Pareto k diagnostic values:
                         Count Pct.    Min. ESS
(-Inf, 0.7]   (good)     5513  98.7%   876     
   (0.7, 1]   (bad)        56   1.0%   <NA>    
   (1, Inf)   (very bad)   14   0.3%   <NA>    
See help('pareto-k-diagnostic') for details.

Code

plot(loo_PH_withC)

if (any(pareto_k_values(loo_PH_withC) > 0.7)) {
  loo_PH_withC_clean <- loo(log_likelihood_PH_withC, k_threshold = 0.7)
  print(loo_PH_withC_clean)
  plot(loo_PH_withC_clean)
}

Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

Computed from 36000 by 5583 log-likelihood matrix.

         Estimate   SE
elpd_loo  -1377.0 38.5
p_loo       152.2 10.9
looic      2754.0 77.0
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume independent draws (r_eff=1).

Pareto k diagnostic values:
                         Count Pct.    Min. ESS
(-Inf, 0.7]   (good)     5513  98.7%   876     
   (0.7, 1]   (bad)        56   1.0%   <NA>    
   (1, Inf)   (very bad)   14   0.3%   <NA>    
See help('pareto-k-diagnostic') for details.

Code

# loo_PH_withC_DU <- loo(log_likelihood_PH_withC_DU)
# print(loo_PH_withC_DU)
# plot(loo_PH_withC_DU)
# 
# if (any(pareto_k_values(loo_PH_withC_DU) > 0.7)) {
#   loo_PH_withC_DU_clean <- loo(log_likelihood_PH_withC_DU, k_threshold = 0.7)
#   print(loo_PH_withC_DU_clean)
#   plot(loo_PH_withC_DU_clean)
# }

2.4 Model comparison

Code

#comp <- loo_compare(loo_main_clean, loo_main_DU_clean, loo_PH_withC_clean, loo_PH_withC_DU_clean)
comp <- loo_compare(loo_main_clean, loo_main_DU_clean, loo_PH_withC_clean)
print(comp, digits = 3)

       elpd_diff se_diff 
model3    0.000     0.000
model1 -120.218    16.923
model2 -138.104    16.699

2.5 High pareto k-values: Try 10-fold CV instead

Code

# Load log_pd_kfold for different models
log_pd_main <- readRDS(file.path(out_path, "log_pd_kfold_bandit2arm_delta_main_n58.rds"))
log_pd_main_DU <- readRDS(file.path(out_path, "log_pd_kfold_bandit2arm_delta_main_DU_n58.rds"))
log_pd_PH_withC <- readRDS(file.path(out_path, "log_pd_kfold_bandit2arm_delta_PH_withC_n58.rds"))
#env_PH_withC_DU <- new.env()
#log_pd_PH_withC_DU <- load(file.path(out_path, "log_pd_kfold_bandit2arm_delta_PH_withC_DU_n58.RData"),env_PH_withC_DU)[1]

# exclude 1st trial per subject as likelihood is identical and hinders loo estimation (causes Inf pareto k values)
log_pd_main <- log_pd_main[,x] 
log_pd_main_DU <- log_pd_main_DU[,x]
log_pd_PH_withC <- log_pd_PH_withC[,x] 
  
# exclude missing trials
log_pd_main <- log_pd_main[,log_pd_main[1,]!=-999]
log_pd_main_DU <- log_pd_main_DU[,log_pd_main_DU[1,]!=-999]
log_pd_PH_withC <- log_pd_PH_withC[,log_pd_PH_withC[1,]!=-999]
#log_pd_PH_withC_DU <- log_pd_PH_withC_DU[,log_pd_PH_withC_DU[1,]!=-999]

2.5.1 Fit per model

Code

elpd_kfold_main <- elpd(log_pd_main)
elpd_kfold_main_DU <- elpd(log_pd_main_DU)
elpd_kfold_PH_withC <- elpd(log_pd_PH_withC)
#elpd_kfold_PH_withC_DU <- elpd(elog_pd_PH_withC_DU)

print(elpd_kfold_main)

Computed from 36000 by 5583 log-likelihood matrix using the generic elpd function

     Estimate    SE
elpd  -2846.5 106.9
ic     5692.9 213.8

Code

print(elpd_kfold_main_DU)

Computed from 36000 by 5583 log-likelihood matrix using the generic elpd function

     Estimate    SE
elpd  -2757.1 102.8
ic     5514.2 205.5

Code

print(elpd_kfold_PH_withC)

Computed from 36000 by 5583 log-likelihood matrix using the generic elpd function

     Estimate    SE
elpd  -2491.3  82.0
ic     4982.6 163.9

Code

#print(elpd_kfold_PH_withC_DU)

2.5.2 Model comparison

Code

#comp2 <- loo_compare(env_main[[elpd_main]], env_main_DU[[elpd_main_DU]], env_PH_withC[[elpd_PH_withC]], env_PH_withC[[elpd_PH_withC_DU]])
comp2 <- loo_compare(elpd_kfold_main, elpd_kfold_main_DU, elpd_kfold_PH_withC)
print(comp2, digits = 3)

       elpd_diff se_diff 
model3    0.000     0.000
model2 -265.796    47.870
model1 -355.182    48.937

3 Step 2: Compare PH withC with different hierarchical structures

3.1 Load Stan fit file

Code

fit_file_PH_withC_hierarchical_group <- 'fit_n58_2024-04-12_bandit2arm_delta_PH_withC_hierarchical_group_estimation1_delta0.9_stepsize0.5.rds'

fit_PH_withC_hierarchical_group <- readRDS(file.path(out_path, fit_file_PH_withC_hierarchical_group))

3.2 Extract log likelihood

Code

# extract log likelihood for each choice
log_likelihood_PH_withC_hierarchical_group <- extract_log_lik(fit_PH_withC_hierarchical_group, parameter_name = "log_lik", merge_chains = T)

# create logical vector coding if a column should be included in log_likelihood
x <- logical(stan_data_hierarchical$T) # create vector containing total number of trials * FALSE
a <- c(stan_data_hierarchical$N+1:(stan_data_hierarchical$N*49)) 
b <- c(stan_data_hierarchical$N*50+stan_data_hierarchical$N+1:(stan_data_hierarchical$N*49))
c <- c(a,b) # concatenate both number vectors
x[c] <- TRUE # set logical value to TRUE if index in c

# exclude 1st trial per subject as likelihood is identical and hinders loo estimation (causes Inf pareto k values)
log_likelihood_PH_withC_hierarchical_group <- log_likelihood_PH_withC_hierarchical_group[,x] 

# exclude missing trials
log_likelihood_PH_withC_hierarchical_group <- log_likelihood_PH_withC_hierarchical_group[,log_likelihood_PH_withC_hierarchical_group[1,]!=-999]

3.3 Fit per model

Code

loo_PH_withC_hierarchical_group <- loo(log_likelihood_PH_withC_hierarchical_group)

Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

Code

print(loo_PH_withC_hierarchical_group)

Computed from 36000 by 5583 log-likelihood matrix.

         Estimate   SE
elpd_loo  -1428.3 40.8
p_loo       117.5  9.3
looic      2856.7 81.7
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume independent draws (r_eff=1).

Pareto k diagnostic values:
                         Count Pct.    Min. ESS
(-Inf, 0.7]   (good)     5532  99.1%   946     
   (0.7, 1]   (bad)        41   0.7%   <NA>    
   (1, Inf)   (very bad)   10   0.2%   <NA>    
See help('pareto-k-diagnostic') for details.

Code

plot(loo_PH_withC)

if (any(pareto_k_values(loo_PH_withC_hierarchical_group) > 0.7)) {
  loo_PH_withC_hierarchical_group_clean <- loo(log_likelihood_PH_withC_hierarchical_group, k_threshold = 0.7)
  print(loo_PH_withC_hierarchical_group_clean)
  plot(loo_PH_withC_clean)
}

Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

Computed from 36000 by 5583 log-likelihood matrix.

         Estimate   SE
elpd_loo  -1428.3 40.8
p_loo       117.5  9.3
looic      2856.7 81.7
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume independent draws (r_eff=1).

Pareto k diagnostic values:
                         Count Pct.    Min. ESS
(-Inf, 0.7]   (good)     5532  99.1%   946     
   (0.7, 1]   (bad)        41   0.7%   <NA>    
   (1, Inf)   (very bad)   10   0.2%   <NA>    
See help('pareto-k-diagnostic') for details.

3.4 Model comparison

Code

comp3 <- loo_compare(loo_PH_withC_clean, loo_PH_withC_hierarchical_group_clean)
print(comp3, digits = 3)

       elpd_diff se_diff
model1   0.000     0.000
model2 -51.328    56.155

3.5 High pareto k-values: Try 10-fold CV instead

Code

# Load log_pd for model with fixed effect for group
log_pd_PH_withC_hierarchical_group <- readRDS(file.path(out_path, "log_pd_kfold_bandit2arm_delta_PH_withC_hierarchical_group_n58.rds"))

# exclude 1st trial per subject as likelihood is identical and hinders loo estimation (causes Inf pareto k values)
log_pd_PH_withC_hierarchical_group <- log_pd_PH_withC_hierarchical_group[,x] 

# exclude missing trials
log_pd_PH_withC_hierarchical_group <- log_pd_PH_withC_hierarchical_group[,log_pd_PH_withC_hierarchical_group[1,]!=-999]

3.5.1 Fit per model

Code

elpd_kfold_PH_withC_hierarchical_group <- elpd(log_pd_PH_withC_hierarchical_group)

3.5.2 Model comparison

Code

comp4 <- loo_compare(elpd_kfold_PH_withC, elpd_kfold_PH_withC_hierarchical_group)
print(comp4, digits = 3)

       elpd_diff se_diff 
model2    0.000     0.000
model1 -114.177   112.101